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We develop a scheme for analytic continuation of the strong-coupling perturbation series of the 
pure Bose-Hubbard model beyond the Mott insulator-to-superfluid transition at zero temperature, 
based on hypergeometric functions and their generalizations. We then apply this scheme for com¬ 
puting the critical exponent of the order parameter of this quantum phase transition for the two- 
dimensional case, which falls into the universality class of the three-dimensional XY model. This 
leads to a nontrivial test of the universality hypothesis. 
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I. INTRODUCTION 


In many branches of theoretical physics one encoun¬ 
ters the necessity to reconstruct an observable from a 
diverging perturbation series. This is possible because 
a divergent series, by itself, is not a bad thing, but still 
carries profound mathematical meaning [1, 2j. Take, for 
example, the geometric series 


designing rapidly converging strong-coupling expansions 
from divergent weak-coupling perturbation series HE]- 
When resorting instead to the more familiar Pade ap¬ 
proximation technique mi, one introduces rational ap- 
proximants of the form 


Al/m(X) 


Em=0 Pn^ n 

i + Ef=i^ 


(3) 



oo 


X- 

n—0 


(i) 


where the sum on the right-hand side (r.h.s.) converges 
only for |z| < 1. Nonetheless, if one regards the left- 
hand side (l.h.s.) as a definition of the formal sum, 
as already done by Euler pQ, that sum obtains a well- 
defined meaning even for \z\ > 1, giving, for instance, 
1 + 2 + 4 + 8 + ... = —1. Similarly, quantum-mechanical 
perturbation theory may yield a formal series 


A(\) ~ ]T a„A" (2) 

n=0 


for some quantity H(A), where A is a small parame¬ 
ter. The coefficients a n may be such that the series is 
asymptotic, having zero radius of convergence, as it hap¬ 
pens, e.g., when computing the ground-state energy of 
an anharmonic oscillator with a quartic perturbation of 
a quadratic potential [3]. The task then again is to iden¬ 
tify the unknown true observable on the l.h.s. from the 
given formal sum on the r.h.s. 

The concept that naturally comes into play here is 
analytic continuation. Still, for applying this concept 
in practice, when merely a few leading coefficients a n 
are available, one needs to invoke some sort of a pri¬ 
ori hypothesis about H(A), either explicitly or implic¬ 
itly. For instance, if one possesses explicit knowledge 
of .A(A) for large values of A, one can exploit this for 
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and equates the coefficients obtained from a Taylor series 
expansion of A L / M up to the order 0(X M+N ) to those of 
the perturbation series (2]). While the resulting Pade ta¬ 
ble then may yield good numerical values of the desired 
quantity, one is implictly imposing the asymptotic be¬ 
havior A l / m { A) ~ pl\ l ~ m /qm for large A, which may 
not be physically correct. 

Recently, Mera, Pedersen, and Nikolic [5] have sug¬ 
gested to replace the rational Pade approximants ^ 
by hypergeometric functions, which represent another 
form of an implicit a priori hypothesis. The examples 
from elementary single-particle quantum mechanics stud¬ 
ied by these authors suggest that the corresponding ana¬ 
lytic continuation technique can dramatically outperform 
Pade and Borel-Pade approaches. Hence, it was con¬ 
jectured that the hypergeometric-function scheme might 
also be useful for many-body problems of condensed- 
matter physics [SJ. 

In the present letter we provide first evidence which 
strongly supports this conjecture. We consider the two- 
dimensional (2d) Bose-Hubbard model on a square lat¬ 
tice, which constitutes a system of paradigmatic impor¬ 
tance for quantum many-body theory maun, and de¬ 
velop a “hypergeometric” technique for obtaining the 
critical exponents of its Mott insulator-to-superfluid 
transition. 

We proceed as follows: We first briefly introduce the 
Bose-Hubbard model, and the strong-coupling perturba¬ 
tion series derived from it, which serves as input for the 
subsequent analysis. Next, we show how Gaussian hyper¬ 
geometric functions 2 -P 1 , and their generalizations q +\F q , 
emerge quite naturally when studying the quantum phase 
transition. We then apply our scheme for computing the 
phase diagram and the critical exponent of the order pa- 
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rameter. This is of some conceptual interest, since these 
exponents are supposed to be universal, but there is still 
a slight discrepancy [ T5 l between experimentally mea¬ 
sured jTBj and theoretically calculated El values. Our 
scheme opens up a fresh approach to this subtle issue. 

II. THE MODEL 


where the expectation value is taken with respect to the 
ground state of the extended model ([ 8 ]), and M denotes 
the number of sites, assumed to be so large that finite-size 
effects do not matter. From this we obtain the suscepti¬ 
bility 

B£ 

2 ^ = — = 2 ( 6 ,), ( 10 ) 


We consider the pure Bose-Hubbard model for inter¬ 
acting Bose particles EMU] on a two-dimensional square 
lattice at zero temperature. In its grand-canonical ver¬ 
sion this model is characterized by three parameters: The 
hopping matrix element J, which quantifies the strength 
of the tunneling contact between neighboring lattice sites, 
the repulsion energy U provided by each pair of particles 
occupying the same lattice site, and the chemical poten¬ 
tial /i. For a given value of /x the competition between 
particle delocalization due to tunneling and localization 
caused by repulsion leads to the well-known quantum 
phase transition from a Mott insulator to a superfluid 
when the ratio J/U is gradually increased, starting from 
zero [10]. We employ the Fock-space operators b\ and b i 
which create or annihilate a Boson at the xth site, so that 


where the first identity constitutes the definition of tp, 
and the second is provided by the Hellnrann-Feynman 
theorem. When taken at rj = 0, this derivative describes 
the response of the original Bose-Hubbard model ([5]) to 
the sources and drains: The expectation value ( 6 ,) = ip is 
zero in the Mott-insulating phase, but takes on nonzero 
values in the superfluid phase, and thus serves as order 
parameter. 

Assuming now that the ground-state energy per site 
can be expanded in a power series of 77 , we write 

£(J/U,n/U,rj) 

oo 

= e 0 {J/U, n/U) + Y, c 2fc (■ J/U, n/U) v 2k ■ (11) 

k =1 


= b V b i ( 4 ) 

counts the number of particles at that site, and use U as 
the energy scale of reference. The non-dimensionalized 
Hamiltonian then is written as 


Hbh = H 0 + H tun , (5) 

where the site-diagonal part 

#o = ( 6 ) 

i i 

models the on-site repulsion and incorporates the chem¬ 
ical potential to fix the particle number; this opera¬ 
tor ([ 6 ]) serves as the starting point for perturbative ex¬ 
pansions HU The further term 

H tnn = -J/uY b \ b j ( 7 ) 

(hi) 

accounts for the tunneling effect, with the sum ranging 
over all pairs of neighboring sites i and j. Following 
Ref. [13j . we then break the particle-number conserva¬ 
tion implied by this Bose-Hubbard model ([5]) by adding 
spatially uniform sources and drains, 

H = H BH + Yv( b \+ b i), ( 8 ) 

i 

where, without loss of generality, we have taken the di¬ 
mensionless source strength 77 to be real. The quantity 
of interest now is the intensive ground-state energy 

£(j/U,fi/U,r,) = (H)/M , (9) 


For each \i/U the coefficients c 2 fc, known as ^-particle 
correlation functions, are then expanded in powers of 

J/U: 

OO 

c 2fe (J/U, ijl/U) = Y (m/£7) {J/U) V ■ (12) 

v=0 

In order to make contact with the Landau theory of phase 
transitions [T5H2(J] . the key idea now is to employ ip in¬ 
stead of 77 as independent variable. This is achieved by 
means of a Legendre transformation, which leads to the 
effective potential hheo 


F = £- 2 ipr) 

= e 0 + ci 2 ip 2 + d^ip 4, + aeip 6 + 0(ip 8 ) (13) 


with the one-particle-irreducible vertices 


1 C4 

a2 =-) «4 = -4 ) «6 

c 2 C 2 



(14) 


having suppressed their dependence on J/U and n/U. 
Since rj and ip constitute a Legendre-conjugated pair, this 
construction implies 


5T 

84’ 


= 2?7 , 


(15) 


leading to the physical interpretation of the formalism: 
Since the actual Bose-Hubbard system ([ 5 ]) is recovered 
by setting 77 = 0 , the physical solutions correspond to 
the stable stationary points of T Ha is]. 

Now one can invoke a standard argument: Assuming 
a 4 and ae to be positive, and neglecting higher order 
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terms of the expansion (13), a single minimum of T is 


found at = 0 as long as 02 > 0, indicating the Mott 
insulator phase. In contrast, if 02 < 0 the minimum is 
found at ip m i n yf 0 , thus signaling the presence of the 
superfluid phase. Therefore, for given chemical potential 
/i/C/ the transition occurs when 02 = — I/C 2 = 0 , that is, 
at that value ( J/U) c at which the series 


C 2 (j/U,fl/U) =Y / ^2^/U) ( J/U) U 


(16) 


,=0 


starts to diverge 
dau form T ss eg + a 2 ip 2 


Moreover, from the usual Lan- 
I- one obtains 


V- 2 . ==^ 

2 a 4 


(17) 


for J/U > ( J/U) c ■ Assuming 04 to be positive and 
smooth at the transition, the exponent (3 which charac¬ 
terizes the emergence of the order parameter according 
to 

^Ln^[J/U-(J/U) c }^ (18) 

is thus solely determined by 02 = — l/c 2 - This sets 
the stage for the present work: Its starting point is 
the perturbation series (16) for the coefficient C 2 . Al¬ 


though this series requires a small parameter J/U it is 
referred to as a strong-coupling expansion m , since it 
should converge in the strongly correlated Mott regime. 
We have evaluated its coefficients a 'P numerically up 
to the order ax = 10 in J/U making use of 

the process-chain approach as devised in general form by 
Eckardt [25] . This technique, which has been recognized 
as an extremely powerful method [251 . is based on Kato’s 
non-recursive formulation of the Rayleigh-Schrodinger 
perturbation series m- Here we take these coefficients 
as input for determining optimal hypergeometric approx- 
irnants to the Landau parameter 02 , as detailed in the 
following section, from which the respective exponents /3 
can then be read off directly. 


III. THE METHOD 


Given the coefficients cq'(m/ 17) for v = 0, 1, 2, ..., 
r m ax, the first task is to deduce the radius of convergence 
of the series (16). This can be accomplished only if some 


a priori knowledge concerning the unknown higher-order 
coefficients is invested. To this end, useful guidance is 
provided by the case of high dimensionality d: As ex¬ 
plained in Ref. [23], for d —> 00 the expansion (16) be¬ 
comes a geometric series, 


(0) 

C2 = «2 2^ 

,=0 


(-2 da^ (J/UY , (19) 


from which one can immediately read off its radius of 
convergence 

-1 


(J/U) c = 


2 d a 


(o) 


( 20 ) 


24 


Is 20 


16.74 

16 


ratios cx//^ /cx// 1J of subsequent coefficients 
- limit expected from QMC data 


perturbation order v 


10 


FIG. 1: Ratios a^V a)" ^ of subsequent coefficients of the 
series (161 for d = 2 and fi/U = 0.3769, as corresponding to 
the tip of the lowest Mott lobe shown in Figj3| The solid hor¬ 
izontal line indicates the limit 1 /(J/[/)]? M == 1/0.05974(3) 

expected from QMC calculations m Observe that all co¬ 
efficients a(2 have the same sign, which impedes the Borel 
summability of the series (161. 


after working out a-p . this leads precisely to the mean- 
field phase boundary [Tp, (23] 

» ^ la Z'* 7 la) ' (21) 

where the integer filling factor g satisfies p/U + 1 > g > 
ji/U. Consequently, in this limiting case the Landau 
coefficient 02 = — I/C 2 takes the simple form 


02 = WG 1 1 “ 


J/U 

(jju). 


( 22 ) 


exhibiting the mean-field exponent 2/3 = 1. 

For finite dimension d , however, the ratio a//^/a// ^ 
of subsequent coefficients is not constant. For d = 2 this 
is shown exemplarily in Fig. [l] for \x/U = 0.3769, corre¬ 
sponding to the tip of the lowest Mott lobe of the phase 
diagram depicted later in Fig. [3] Thus, the corresponding 
series (16) is no longer geometric, so that it is tempting 


to assume a Landau coefficient of the form 


a 2 = ( 1 - 


a 


(o) 


J/U 

VJu\ 


2/3 


(23) 


thereby admitting nontrivial exponents 2/3. According to 
this educated guess, the expansion (16) should have the 
form of a binomial series, 


(o) \ ' 

C 2=«2 L 

i *=0 


(2/3), ( J/U 
v\ {(J/U) ( 


(24) 


where (a), = a(a + l) • • • (a+u — 1) is the usual Pochham- 
mer symbol [28] . From an optimal fit of the given coef¬ 
ficients a ^ to this hypothesis one could then determine 
approximate values of the two parameters 2/3 and (J/U) c . 
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But this guess can still be improved: Realizing that 
the binomial series coincides with the function iFo, em¬ 
ploying the nomenclature used for generalized hyperge¬ 
ometric functions [291 . one may generalize the a priori 
ansatz (24) further and require 



where now 2 -Fi(a, b;c;z) denotes the well-known Gaus¬ 
sian hypergeometric function [251 ES], giving us four de¬ 
grees of freedom for a least-square fit to the i/ max pertur¬ 
bative data. The strength of its singularity at the point 
of divergence is given by 

2 ^i (a, fr; c; z) oc _ ^ a+b _ e for ^ ^ 1 , (26) 

from which one finds the exponents 

2/3 = a + b-c . (27) 


Going still one step further, one may replace 2 F 1 by 
a generalized hypergeometric function q +iF q providing 
2 q + 2 degrees of freedom, requiring the evaluation of the 
perturbation series at least to the corresponding order. 
This possibility to perform analytic continuation of a per¬ 
turbation series by means of an analytic function which 
itself is a member of a general familiy of “higher order” 
functions is the core of the proposal made in Ref. ;5]. 
The hypergeometric functions, just as the binomial se¬ 
ries, can adopt complex values beyond their radius of 
convergence. In view of its physical meaning we require 
C 2 , as well as 0 , 2 , to be a real quantity, so that we take 
the real part of q +iF q . Technically, this is achieved by 
computing lim e ^. 0 ( q+ iF q (x + ie) + q+ iF q (x - ie))/2. 

Thus, from the expectation of a nontrivial exponent 
we infer that C 2 has to have an essential singularity at 
(J/U ) c with a strength determining the respective expo¬ 
nent. This is why a +\ F a are suitable approximants: As 
exemplified by Eq. (26), the strength of their singularities 
can be tuned by adjusting their parameters. 


IV. RESULTS 


We have applied this strategy to the series ( [l 6 | for the 
two-dimensional Bose-Hubbard model with 0 < n/U < 
4, having at disposal its coefficients up to ^ max = 
10 l22H2ij . Figure [ 2 ] again displays the ratios / 01 ^ 

of subsequent coefficients for \x/U = 0.3769, now plotted 
vs. the reciprocal order 1/v. In addition, we also indicate 
by continuous lines the ratios resulting from the best fit 


(Q 

a\ = a 


(o) 


( 1 - 405 ), 


16.51" 


(28) 



FIG. 2: Ratios /at^ ' 1 of subsequent coefficients of the 

series (161 for d = 2 and p,/U = 0.3769, plotted vs. the 
reciprocal order 1 /v, together with the corresponding ra¬ 
tios obtained from the fit ( |28| ) to the binomial series, and 
from the fit ( |29| ) to the Gaussian hypergeometric function. 
Again, the horizontal line indicates the limit 1 /(J/U)^ MC = 
1/0.05974(3) expected from QMC calculations [12] . 


to the binomial series hypothesis (24), and from the best 
fit 


a 


(G 

2 



(l.398) ;y (- 0.7684)^ 
v\ ( - 0.7606)^ 


16.61" 


(29) 


to the Gaussian hypergeometric hypothesis (25). Evi¬ 
dently, the quality of these fits is excellent. This al¬ 
lows us to perform reliable extrapolations to v —> 00, 
where the ratios should approach the expected value 
1 /(J/U)Q MG = 1/0.05974(3) known from quantum 
Monte Carlo (QMC) simulations [T2]. In addition, we 
have also fitted the exact coefficients to those of general¬ 
ized hypergeometric functions 3 F 2 and 4 F 3 . 


Performing these procedures for all values of the chem¬ 
ical potential that are of interest, and reading off the re¬ 
spective values of (J/U) c , we obtain the system’s phase 
diagram. In Fig. [3] we show this diagram for 0 < /i/I/ < 
4, as resulting from the binomial and from the Gaussian 
hypergeometric fit, respectively. The relative deviation 
between both curves stays below 4%. It is largest halfway 
between the position of a tip of a Mott lobe and the near¬ 
est integer values of fi/U\ at the tips these deviations are 
smaller than 1 %. 


The tips of the Mott lobes represent multicriti- 
cal points with particle-hole symmetry; here the sys¬ 
tem falls into the universality class of the (d + 1 )- 
dimensional XY model m- The critical scaled hopping 
strength ( J/U) c at the tip of the lowest lobe for d = 2, 
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- binomial series 1 Fq extrapolation 

-hypergeometric function 2 F\ extrapolation 



FIG. 3: Zero-temperature phase diagram of the 2d Bose- 
Hubbard model as obtained from the binomial series hypothe¬ 
sis, and from the Gaussian hypergeometric hypothesis. Inside 
the lobes the system is in a Mott insulator state with g par¬ 
ticles per site, outside these lobes in a superfluid state. Also 
shown is the expression for the lowest Mott lobe stated in 
Ref. HH. 


as provided by the respective fit, figures as 


iF 0 

(J/U) c 

= 0.06056 

2 -Fl 

(J/U) c 

= 0.06021 

3^2 

(J/U) c 

= 0.06003 

4 F 3 

(J/U) c 

= 0.06004, 


(30) 


which matches the QMC value (J/U)^ MG = 0.05974(3) 
quite well. If we assume this value to be exact, and take 
the result provided by 2 -F 1 as sound compromise between 
the number of coefficients available and the number of fit 
parameters, the error of that result is less than 1%. 

So far, we have used hypergeometric continuation 
merely to reproduce existing knowledge, thus confirm¬ 
ing its reliability. However, with the determinantion of 
the order parameter’s exponent we now move to a ground 
which is technically far more demanding [30] . In Fig. [4] 
we show the Gaussian hypergeometric description of the 
Landau coefficient < 22 , and its analytic continuation be¬ 
yond the transition point, again at the tip of the lowest 
Mott lobe. The data are well described by the power-law 
fit 

8.679((J/[/) c — J/[/) 1 ' 39 °for J/U < ( J/U) c 

/ \ 1 SQD 

-2.945 (J/U - (J/U ) c ) for J/U > ( J/U) c . 

Performing this procedure within the interval 0 < 
n/U < 4.0 for both the binomial and the Gaussian hy¬ 
pergeometric ansatz, we obtain the exponents displayed 
in Fig. [5} While both variants lead to notably different 



FIG. 4: Behavior of the Landau coefficient 02 for [i/U = 
0.3769, as obtained by Gaussian hypergeometric continuation; 
for J/U > (J/U )o the real part is shown. Also shown is the 
power-law fit (|31|). 



exponent of the order parameter ip 


FIG. 5: Exponents /3 of the order parameter ip for the Mott 
insulator-to-superfluid transition, as obtained by binomial 
and Gaussian hypergeometric continuation, respectively. The 
“critical” value « 0.69 is adopted at the tips of the Mott 
lobes. 


results for general n/U, they agree quite well at the tip 
of the lobes, i.e., at the critical points. This observation 
is quite significant, since one expects nontrivial critical 
exponents only at the tips of the lobes, while the system 
should be mean-field like for all other \xjU. Our approx¬ 
imate scheme can only yield continuous lines, and it is 
an open question whether these would reduce to <5-like 
spikes with higher orders. However, the relative stability 
of the data at the lobes’ tips indicates that one can actu¬ 
ally determine the true critical exponent /3 C of the Mott 
insulator-to-superfluid transition with good accuracy by 
hypergeometric continuation. In particular, at the tip of 
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the lowest lobe we obtain the following values: 


^0 

/3c 

= 0.7023 

Fi 

/3 C 

= 0.6949 

f 2 

/3 C 

= 0.6881 

f 3 

/3 C 

= 0.6904 


This gives rise to a nontrivial test of the universality 
hypothesis for critical phenomena. Using the scaling 
relation /? c = (1 + rf)v, and inserting the best known 
estimates for the critical exponents p = 0.0380(4) and 
v = 0.67155(27), as derived from the three-dimensional 
XY model by combining Monte Carlo simulations based 
on finite-size scaling methods and lrigh-temperature ex¬ 
pansions m, the assumption of universality yields the 
expectation /3 C = 0.6971(6). Indeed, this coincides within 
less than 1% with our 2 Tj-estimate extracted from the 
2d Bose-Hubbard model. It remains to be seen whether 
further refinement of our approach will result in an even 
better confirmation of the universality hypothesis. 


V. CONCLUSIONS 


In summary, we haven taken up an idea put forward 
by Mera, Pedersen, and Nikolic, who have suggested to 
utilize hypergeometric functions for the analytic contin¬ 
uation of divergent perturbation series [5]; here we have 
adapted this concept to the strong-coupling perturba¬ 
tion series (16) of the Bose-Hubbard model. After eval¬ 


uating this series to the maximum accessible order in 
the scaled hopping strength J/U , which is ^ max = 10 
in the present case, we are in a position to determine 
the parameters of its hypergeometic approximants from 
a least-square fit with high accuracy. This has enabled 
us to assess the critical exponent of the order parame¬ 
ter of the Mott insulator-to-superfluid transition. Com¬ 
pared to a previous attempt to deduce critical exponents 
from diverging perturbation series |24] , the present ap¬ 
proach is conceptually simpler, and more easy to handle 
in practice. The success of this approach indicates that 
hypergeometric functions, and their generalizations, in¬ 
deed embody the proper a priori knowledge required by 
this quantum phase transition. Aside from further refine¬ 
ments, the next steps to be taken with hypergeometric 
analytic continuation will involve the investigation of the 
superfluid density, and the corresponding analysis of the 
3d Bose-Hubbard model. 
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